function estimatedDimensionality = test_dimensionality(data, dimensionality, includeSlowMethods)
  methods = {'CorrDim' 'NearNbDim' 'EigValue' 'MLE'};
  
  slowMethods = {'PackingNumbers' 'GMST'};
  
  if ~exist('dimensionality', 'var')
    % No intrinsic dimensionality has been suggested
    dimensionality = Inf;
  end
  
  if exist('includeSlowMethods', 'var')
    disp('Including slow methods…')
    
    methods = [methods slowMethods];
  end
  
  totalTime = 0;
  totalDimensionality = 0;
  
  results = {};
  
  disp('Estimating dimensionality…')
  
  for i = 1:numel(methods)
    method = methods{i};
    tic();
    
    disp('')
    disp('---')
    disp(['  ** ' method ' **'])
    fflush(stdout);
    
    estimatedDimensionality = intrinsic_dim(data, method);
    totalDimensionality += estimatedDimensionality;
    
    time  = toc();
    if (dimensionality != Inf)
      error = abs(estimatedDimensionality - dimensionality) / dimensionality;
    else
      error = 0;
    end
    
    totalTime += time;
    
    disp(['  Error: ' num2str(error)])
    disp(['  Time:  ' num2str(time) ' s'])
    disp('---')
    
    result.name  = method;
    result.dim   = estimatedDimensionality;
    result.error = error;
    result.time  = time;
    results      = [results result];
  end
  
  disp('')
  
  fastestMethod.name  = 'None';
  fastestMethod.dim   = Inf;
  fastestMethod.error = Inf;
  fastestMethod.time  = Inf;
  
  bestMethod.name     = 'None';
  bestMethod.dim      = Inf;
  bestMethod.error    = Inf;
  bestMethod.time     = Inf;
  
  for i = 1:numel(results)
    result = results(i);
    
    if (result.error < bestMethod.error)
      bestMethod.name  = result.name;
      bestMethod.dim   = result.dim;
      bestMethod.error = result.error;
      bestMethod.time  = result.time;
    end
    
    if (result.time < fastestMethod.time)
      fastestMethod.name  = result.name;
      fastestMethod.dim   = result.dim;
      fastestMethod.error = result.error;
      fastestMethod.time  = result.time;
    end
  end
  
  if (dimensionality != Inf)
    disp('---')
    disp(['** Best Method:    ' bestMethod.name])
    disp(['   Dimensionality: ' num2str(bestMethod.dim)])
    disp(['   Rel Error:      ' num2str(bestMethod.error)])
    disp(['   Time:           ' num2str(bestMethod.time) 's'])
    disp('---')
  
    disp('')
  end
  
  disp('---')
  disp(['** Fastest Method: ' fastestMethod.name])
  disp(['   Dimensionality: ' num2str(fastestMethod.dim)])
  disp(['   Rel Error:      ' num2str(fastestMethod.error)])
  disp(['   Time:           ' num2str(fastestMethod.time) 's'])
  disp('---')
  
  disp('')
  
  estimatedDimensionality = totalDimensionality / numel(methods);
  
  disp(['Dimensionality: ' num2str(round(estimatedDimensionality))])
  disp(['Total time: ' num2str(totalTime) 's'])
  disp('')
  
